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Finite quantum environments as thermostats: an analysis based 
on the Hilbert space average method 
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Abstract. We consider discrete quantum systems coupled to finite environments which may possibly consist 
of only one particle in contrast to the standard baths which usually consist of continua of oscillators, spins, 
etc. We find that such finite environments may, nevertheless, act as thermostats, i.e., equilibrate the system 
though not necessarily in the way predicted by standard open system techniques. Thus, we apply a novel 
technique called the Hilbert space Average Method (HAM) and verify its results numerically 
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1 Introduction 

Due to the linearity of the Schrodinger equation concepts 
like ergodicity or mixing are strictly speaking absent in 
quantum mechanics. Hence the tendency towards equilib- 
rium is not easy to explain. However, except for some ideas 
[H, 0] the approaches to thermalization in the quantum 
domain seem to be centered around the idea of a ther- 
mostat, i.e., some environmental quantum system (bath, 
reservoir) , enforcing equilibrium upon the considered sys- 
tem. Usually it is assumed that the classical analogon of 
this bath contains an infinite number of decoupled degrees 
of freedom. 

Theories addressing such scenarios are the projection 
operator techniques (such as Nakajima-Zwanzig or the 
time-convolutionless method see Q) and the path integral 
technique (Feynman Vernon 1). The projection operator 
techniques are exact if all orders of the system-bath inter- 
action strength are taken into account which is practically 
unfeasible. However, assuming weak interactions and ac- 
cordingly truncating at leading order in the interaction 
strength, which goes by the name of "Born approxima- 
tion" (BA), produces an exponential relaxation behavior 
(cf. 0, Q) whenever the bath consists of an continuum of 
oscillators, spins, etc. The origin of statistical dynamics 
is routinely based on this scheme, if it breaks down no 
exponential thermalization can a priori be expected. 

In contrast to the infinite baths which are extensively 
discussed in the mentioned literature, we will concentrate 
in this article on finite environments. The models we ana- 
lyze may all be characterized by a few-level-system (S, 
considered system) coupled to a many-level-system (E, 
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environment) consisting of several relevant energy bands 
each featuring a number of energy eigenstates (e.g. see 
Fig. [1]). Thus, this may be viewed as, e.g., a spin cou- 
pled to a single molecule, a one particle quantum dot, an 
atom or simply a single harmonic oscillator. Note that 
the spin, unlike in typical oscillator baths or the Jaynes- 
Cummings Model, is not supposed to be in resonance with 
the environments level spacing but with the energy dis- 
tance between the bands. There are two principal differ- 
ences of such a finite environment level scheme from the 
level scheme of, say, a standard oscillator bath: i) The to- 
tal amount of levels within a band may be finite, ii) Even 
more important, from e.g. the ground state of a standard 
bath there are infinitely many resonant transitions to the 
"one-excitation-states" of the bath. But from all those, the 
resonant transitions lead back to only one ground state. 
Thus, the relevant bands of any infinite bath would con- 
sist of only one state in the lowest band and infinitely 
many states in the upper bands. In this paper we focus 
on systems featuring arbitrary numbers of states in any 
band. (For a treatment of finite baths under a different 
perspective, see @, 111)- 

It turns out that for the above mentioned class of 
models standard methods do not converge and thus the 
unjustified application of the BA produces wrong results 
(cf. [I Q3). This holds true even and especially in the 
limit of weak coupling and arbitrarily dense environmen- 
tal spectra. Nevertheless, as the application of the Hilbert 
Space Average Method (HAM) predicts, a statistical re- 
laxation behavior can be induced by finite baths. It sim- 
ply is not the behavior predicted by the BA. Thus, the 
principles of statistical mechanics in some sense apply be- 
low the infinite particle number limit and beyond the BA. 
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Fig. 1. Two- level system coupled to a finite environment. Spe- 
cial case for only two bands in the environment 



This also supports the concept of systems being driven 
towards equilibrium through increasing correlations with 
their environments [1, [ll], Q3, EH, rather than the idea 
of system and environment remaining factorizable, which 
is often attributed to the B A [3|, |4| . 

Our paper is organized as follows: First we introduce 
our class of finite environment models and the appropriate 
variables (Sect. l2~Tl Sect. 12. 2( ). Then we compute the short 
time dynamics of those variables fSect. [2~3|) . Hereafter we 
introduce HAM and show in some detail how it can be 
exploited to infer the typical full time relaxation from the 
short time dynamics (Sect. [3]). The theory is then veri- 
fied by comparing the HAM predictions with numerically 
exact solutions of the time dependent Schrodinger equa- 
tion for the respective models (Sect. In the following 
Section the limits of the applicability of HAM which turn 
out to be the limits of the statistical relaxation itself are 
discussed (Sect. [5]). Finally we conclude (Sect. [6]). 



2 System and Dynamics 
2.1 Finite Environment Model 

As just mentioned we analyze a few-level-system S, with 
state space Hs, coupled to a single-many-level system E 
with state space He consisting of energy bands, featur- 
ing, for simplicity, the same width and equidistant level 
spacing. A simple example is depicted in Fig. [1] with two 
bands in the environment only. The Hilbert space of states 
of the composite system is given by the tensor product 
H = H s ®He- 

In Tig let us introduce standard transition operators 
Pij — K)0'|, where \j) are energy eigenstates of the 
considered system S. Furthermore, we define projection 
operators that implement projections onto the lower, re- 
spectively upper band of the environment in He by 



n a = S,\n a ){n a \ 



(1) 



where \n a ) are energy eigenstates of E, a labels the band 
number. Those projectors meet the standard property 

n a n a > = 5 a a>n a , . (2) 

Thus, the number of eigenstates in band a is given by 
N a = Tr{77 a }. 



The complete Schrodinger picture Hamiltonian of the 
model consists of a local and an interaction part H = 
H\ oc + V, where the local part reads 



#ioc = H s ® 1 E + Is <8> H E 



(3) 



with 



Hs — E EiPa 



Sen a 
N a 



\n a )(n a \ 



(4) 



Here we introduced the energy levels Ei of S and the 
mean band energies E a of the environment. Note that 
[n a ,H^\ = 0. For the special case depicted in Fig. Q] one 
gets a = 1, 2 only. 

The interaction may, in principle, be any Hermitian 
matrix defined on H. We choose to decompose it uniquely 
as follows 

V = Y t P ij ®C ij , (5) 

ij 

where the (3y themselves may be decomposed as 

Cjj = Cjj,gb with Cij t ab = n a Cijflb ■ (6) 

ab 

For our special case, the interaction V and its decompo- 
sition are sketched in Fig. [2j For later reference we define 
the "coupling strengths" \ij, a b as 



\2 



Km 



(7) 



and due to the Hermiticity of the interaction \ij jCt b = \ji,ba 
is a real number. Conditions on those interaction strengths 
are discussed in more detail in Sec. 15 - 1 1 but in general we 
assume them to be weak compared to local energies in S 
and E, i.e., AE from dJ). 

There are two different types of (additive) contribu- 
tions to V: (7-terms that induce transitions inside the sys- 
tem S (featuring i ^ j) as well as terms which do not 
(featuring i — j). Since the first type exchanges energy 
between system and environment it is sometimes referred 
to as "canonical coupling" V can (those terms are shaded 
grey in Fig. [2]). The second type produces some entangle- 
ment thereby causing decoherence (those terms are white 
in Fig. [2]), but does not exchange energy and therefore 
refers to a "microcanonical coupling" V m - lc in the context 
of quantum thermodynamics (cf. [15| ) . We do not specify 
the interaction in more detail here. To keep the following 
theoretical considerations simple we only impose two fur- 
ther conditions on the interaction matrix. First we require 
that the different parts of V as displayed in Fig. [2] are not 
correlated unless they are adjoints of each other. Hence 
we get for the following traces 

Tr{Cy,ai> Cj'i> tb 'a'} ~ N a N b A?- ab 6i,i>6jj>5 ata >5b,b> ■ (8) 
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Fig. 2. Scheme of the interaction matrix for the special model 
Fig. [T] explanation see text. 



Furthermore, demand the traces over individual contribu- 
tions of V to vanish, i.e., 



Tr{C ij<ab } « . 



(9) 



Both constraints are common in the field and definitely 
apply for the used numerical examples. 

For all numerical investigations (see Sect.0]) we are us- 
ing complex Gaussian distributed random matrices with 
zero mean to model the interaction. Thus, the mentioned 
special conditions apply: Only adjoint blocks are corre- 
lated and the above traces are extremely small for Gaus- 
sian random numbers with zero mean. This interaction 
type has been chosen in order to keep the model as gen- 
eral and free from peculiarities as possible. For example, 
in the fields of nuclear physics or quantum chaos random 
matrices are routinely used to model unknown interaction 
potentials. We do, however, analyze the dynamics gener- 
ated by one single interaction, not the average dynamics 
of an Gaussian ensemble of interaction matrices. 



2.2 Reduced Dynamics and Appropriate Variables 

Of course we are mainly interested in the time evolution 
of the system S separately, i.e., we would like to find an 
autonomous, time-local equation for the dynamics of its 
reduced density matrix p. However, it turns out that an 
autonomous description in terms of p is, in general, not 
feasible for finite environments. We thus aim at finding an 
autonomous set of equations for the dynamics of a set of 
variables that contain slightly more information than p, 
such that from the knowledge of this set p may always be 
computed. 

We simply name the set here, and explain in the follow- 
ing the derivations of its dynamics. Consider the following 



operators 



Pij, a — Pij H a 



According to ([2]) we find 

P ■ P , ; , — r) ■,() ,P, , 



(10) 



(11) 



and thus the given operators together with zero form a 
group which we mention here for later reference. Through- 
out this paper we think of the full system as being always 
in a pure state \ij})(4>\- Thus the expectation values of the 
operators (fTU|) may be denoted as 



\Pii,aW) = P i3,a 



(12) 



For the dynamics of those expectation values we are going 
to derive an autonomous set of equations. In terms of those 
variables the reduced density matrix elements pij read 



Pij 



(13) 



This may simply be computed from the definition of den- 
sity matrix of S 



pa = mi) = WT^mmij) 



(14) 



2.3 Short Time Dynamics 



In order to find the full dynamics of the P's defined in 
(|12|) . we start of by computing their short time evolu- 
tions in this Section. To those ends we change from the 
Schrodinger to the interaction (Dirac) picture in which the 
originally constant interaction V from Sect. |2~T1 becomes 
time dependent 



V(t) = e 



_ JH loc t/h 



Ve 



-iH loc t/h 



(15) 



As well-known in the interaction picture the time evolu- 
tion may be written in terms of an propagator D{t, i) 



m+T))=D{ T ,t) \m) 



(16) 



(from here \ip(t)) refers to the interaction picture). The 
propagator D (t, t) may be explicitly noted in terms of an 
Dyson expansion 



D( T ,t) = i 



E 



where 



U j (r,t)=T'[[ 



r n +t 



dr„ V(t„ + 1) 



(17) 



(18) 



with T being the standard time ordering operator. 

The time evolution of the expectation values Pij, a ac- 
cording to (fT2"l) reads 



PijAt + T ) = (V>(* + r)\Pij, a \ip{t + t)) 



(19) 
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which using (|16|) may also be written as 

p^At +t) = (m\PijAt + r)\m) > (20) 

with 

P^, a {t + r) = D t (T,t)P lJ , a D{T,t). (21) 

The above definition allows to write the expectation value 
of Pij,a at time t + r (in the interaction picture) as an 
expectation value of some operator Pij, a (t + t) at time t. 
This particular form is well suited to asses that dynamics 
by HAM as will be explained in the next Section. 

If t is short and the interaction is weak the propagator 
may be approximated by a truncation of the Dyson series 
(fT7|) to leading order which is in this case second order. 
Let the truncated propagator be denoted as T)% (r, t) . This 
truncated propagator can typically be computed, even if 
complete diagonalization is far beyond reach (for a more 
explicit treatment of cf. App. A). Thus the approxi- 
mate form we are going to use in Sect. 13.21 reads 

Pij, a (t + r) = £>1(t, t) P^ a D 2 (t, t) . (22) 

3 Dynamical Hilbert Space Average Method 

3.1 Definition and Calculation of the Hilbert Space 
Average 

The Hilbert space average method (HAM) is in essence a 
technique to produce guesses for the values of quantities 
defined as functions of a wave function if itself 
is not known in full detail, only some features of it. In 
particular it produces a guess for some expectation value 
(ip\S\ip) [cf. ([20]) ] if the only information about is the 
set of expectation values ( - 0|f > r),a|V ; ) — Pij,a mentioned 
below. Such a statement naturally has to be a guess since 
there are in general many different that are in accord 
with the given set of Py, a , but produce possibly different 
values for (ip\S\ip). The question here is whether the distri- 
bution of (4>\S\ip) , s produced by the respective set of Ws 
is broad or whether almost all those \tp)'s yield (ip\S\ip)'s 
that are approximately equal. It turns out that if the spec- 
tral width of S is not too large and S is high-dimensional 
almost all individual \ip) yield an expectation value close 
to the mean of the distribution of (V>|5'|V ; ) ' s (see Sect. [5] 
and [15]). The occurrence of such typical values in high- 
dimensional systems has recently also been exploited to 
explain the origin of statistical behavior in (l6l . [l7| . 

To find the above mean one has to average with respect 
to the |V>)'s. We call this a Hilbert space average S and 
denote it as 

S=msm m p ij ^ Pijia y (23) 

This expression stands for the average of (ip\S\ip) over all 
l^) that feature {^}\Pij ia \ijj) = Pij,a but are uniformly dis- 
tributed otherwise. Uniformly distributed means invari- 
ant with respect to all unitary transformations e that 



leave the respective set of expectation values unchanged, 

i.e., (?/)|e i6 Pij,ae" i6 |-0) = (ip\Pij,a\ip) • Thus the respective 
transformations may be characterized by 

[G, 4>] = . (24) 

Instead of computing the so defined Hilbert space average 
(|2"3")) directly by integration as done, e.g. in [15|, [lj| we will 
proceed in a slightly different way, here. To those ends we 
change from the notion of an expectation value of a state 
to one of a density operator 

s=msm = i^{smm, (25) 

where we skipped the constant expectation values of the 
Hilbert space average for the moment. Exchanging the 
average and the trace, one may rewrite 

S = Tr{S[|V>>M]}=Tr{Sa} (26) 

with 

a = Mm m p ijM=Pi . a} . (27) 

To compute a we now exploit its invariance properties. 
Since the set of all that "make up" a [that belong 
to the averaging region of (|27|) ] is characterized by being 

invariant under the above transformations e _ , a itself 
has to be invariant under those transformations, i.e. 

e i6 ae- id = a. (28) 

This, however, can only be fulfilled if [G, a] = for all 
possible G. Due to (f2~4")l the most general form of a which 
is consistent with the respective invariance properties is 

& = ^2Pij,a.Pij,a , (29) 
ija 

where the coefficients Pij,a are still to be determined. In 
principle the above sum could contain addends of higher 
odcr, i.e., products of the P-operators, but according to 
the properties of the projection and transition operators 
[especially (ITTI) ]. those products reduce to a single P-op- 
erator or zero (in other words, the Pij, a form a group), 
hence (|29|) is indeed the most general form. 

How are the coefficients Pij t a to be determined? From 
the definition of a in (|27|) it follows 

Tr{&P ifj , ia ,} = P vj , ta , . (30) 

By inserting (|29|) into ([30|) and exploiting (|11| the coeffi- 
cients are straightforward found to be 

/',.,• = ■ ( 31 ) 

Thus, wc finally get for the Hilbert space average (f26)) 

S = Tr{5d} = £ %Tr{.S7 ',,,!• . (32) 
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Fig. 3. Repeated guessing scheme. To each iteration time step 
t corresponds an increasing uncertainty (variance) of the guess. 



For P's featuring i = j one finds 
Pu,a(* + T)-P iij0 (t) 

Pmm,b{&) Pit, a (^) 



^ 2Re/ im , g& (T) ( - 

mb 



Nh 



and for the P's with i =/= j 
P ij<a (t + T)-P ijta (t) 



2 N n 



where the /(r)'s are defined as 



(35) 



(fj?n,ab( T ) + fun,ab( T )) i (36) 



3.2 Iterative Guessing 

To find the (reduced) autonomous dynamics for the Pji, a 
from HAM we employ the following scheme: Based on 
HAM we compute a guess for the most likely value of the 
set Pji t a at time (t + r) [i.e.,P,-j )0 (t + r)] assuming that we 
knew the values for the Pji, a at time t [i.e.,Pji >a {t)]. Once 
such a map Pji, a (t) — > Pji, a (t + t) is established it can 
of course be iterated to produce the full time dynamics. 
This of course implies repeated guessing, since in each iter- 
ation step the guess from the step before has to be taken 
for granted. However, if each single guess is sufficiently 
reliable, i.e., the spectrum of possible outcomes is rather 
sharply concentrated around the most frequent one (which 
one guesses), even repeated guessing may yield a good 
"total" guess for the full time evolution. The scheme is 
schematically sketched in Fig. [3] Some information about 
the reliability of HAM guesses has already been given in 
Sect. 13.11 the applicability of the whole scheme will be 
analyzed more thoroughly in Sect. [51 

To implement the above scheme we consider the equa- 
tion one gets from inserting Py, a (i + r) for S in (|32|) 



Pij,a( t + T 



WlA'/,a'W=iV4',a'(*)} 



i'j'a' 



P °' l ' N a ' i (t) Tr{P^(r, t)P ijja D 2 (T, t)A'i',,'} . (33) 



This is the Hilbert space average (HA) over all possible 
Pji,a(t + r) under the condition that one had at time t the 
set Pji, a {t). Thus the (iterative) guess now simply consists 
of replacing the HA by the actual value, i.e., 



i'j'a' 



The evaluation of the right hand side requires some rather 
lengthy calculations, but can be done without further as- 
sumptions or approximations. The interested reader may 
find the details in App. B. Here we simply give the results 
and proceed. 



ikafc(r) := / dr' / Ar" 9ij , ab {r") (37) 
Jo Jo 

<fea 6 (r") := ^ Tr E {C aMj (T") C ba<ji } . (38) 

Note that ([35|) and (|36f now are autonomous (closed) in 
terms of the respective P's and there is no more explicit 
dependence on the absolute time t. It turns out (see below) 
that the /(t)'s are approximately linear in r. Hence for 
the squared absolute values of the P's with i ^ j (which 
we will be primarily analyzing rather than the P's them- 
selves) one finds to linear order in r 

\Pij,a(t + T)\ 2 -\PijAt)\ 2 

= - feWl? (Re/ J m, ab (r) + Re/ im , a6 (r)) . (39) 



3.3 Correlation Functions and Transition Rates 

In order to interpret (|35[) and (|39[) appropriately, we need 
some information about the correlation functions /(r). 
Apparently those /(t)'s are essentially integrals over the 
same environmental temporal correlation functions g(r") 
that appear in the memory kernels of standard projection 
operator techniques. (Only here they explicitly correspond 
to transitions between different energy subspaces of the 
environment.) Thus we analyze the g(r")'s from (|37|) more 
thoroughly. Their real parts (which eventually essentially 
matter) read 



Regij,ab{T") =ra £ \{n a \Cij\n b )\ 2 x 



(40) 



cos ( \Ei — Ej + E(n a ) — E(rib) 



Thus they simply consist of a sum of weighted cosine func- 
tions with different frequencies. The set of those weights 
essentially gives the Fourier transform of the correspond- 
ing correlation function. First of all, only if the transition 
within the system (j — > i) is in resonance with the en- 
ergy gap between the bands a, b, g(r) will contain any 
small frequency contributions at all. Hence, only in this 
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case temporal integrations, i.e., the corresponding /'s will 
be nonzero. In the resonant case the frequency spectrum 
will stretch from zero to a frequency on the order of 8e/K, 
at least if the interaction gives rise to the corresponding 
transitions of the environment. Thus g(r) will decay on a 
timescale on the order of r r with 



h 

Te 



(41) 



For t > t c , g(r) will be essentially zero. This means that 
/(r) which is a twofold temporal integration of g(r) will 
grow linear in time, i.e., /(r) = jt after t t c . The factor 
7 is given by the area under the curve g(r) up to approx- 
imately t c . If t c was infinite 7 would only be determined 
by the weight of the zero-frequency terms of g{r). Since r c 
is finite, 7 is related to the "peak density" of g within fre- 
quency range from zero to Auj with Auj c <C 1/t c ~ 5e/H. 
Which means 7 is eventually given by the sum of all 
weights that correspond to frequencies from zero to Auj 
divided by Auj and multiplied by it. Since in our model 
the |(n a |Cij are Gaussian distributed random num- 
bers we eventually find for the /(t)'s 



Re/^, a fc(T) 



hSe 



(42) 



This result can apparently be connected to the transition 
rate as obtained from Fermi's Golden Rule. Let jij^b be 
the Golden Rule transition rate for a transition of full 
system characterized by j — > i and b — > a. Then the con- 
nection reads 



2 Re fim,ab(.T~) ~ 'Yim^ab 

N b T. 



(43) 



Since the Golden Rule transition rates depend on the state 
densities around the final states, respective "forward" and 
"backward" rates are, for equal bandwidths, connected as 



Na 



(44) 



3.4 Reduced Equations of Motion 

Inserting (|4*3"|) and (HH) into (|3"5|) and ([55]) allows for a com- 
putation of the full dynamics of the Pu ya and the |-FV/, Q | 2 
through iteration. The iteration has to proceed in time- 
steps that are longer than r c , but shorter than r^. The 
latter will be explained in Sect. 15.21 Assuming that t c is 
short compared to the timescale of the relaxation dynam- 
ics it may be written as 

^rPiiA 1 ) = E [7im,abPmm,b( t ) ~ lmi,ba P H,a{t)] , (45) 



dt 

d_ 

dt 



mb 



PijAVf = -|^,a W| 2 ^2hmi,ba + Imjjba) • ( 46 ) 
mb 



(This form is in accord with recent results from novel pro- 
jection operator techniques Q) We now analyze this set 



of equations in a little more detail. The Pu ya may be inter- 
preted as the probability to find the joint system in state 
i for S and in band a with respect to the environment. 
Equation (|45|) obviously has the form of a master equa- 
tion, i.e, the overall probability is conserved and there is 
a stable fixpoint which sets the equilibrium values for the 
Pu : a- According to (|46|) the Pij. a will all decay to zero. 
Taking (|13p into account this implies that p will reach an 
equilibrium state which is diagonal in the basis of the en- 
ergy eigenstates of S. As already mentioned below (|40| 
transitions occur only between resonant states, i.e., the 
Jim.ab are zero unless E rn + E b « Ei + E a where E a , E b 
are the corresponding mean band energies. Thus, if we 
define the approximate full energy of some state E 



E = En 



E h 



E = Ei + E n 



(47) 



we may label full system states by i, E(m, E) rather than 
i,a(m,b) and nonzero transition rates by im,E rather 
than im,ab, i.e, P iiya — > Pf , Jim,ab —> lf m - With this 
index transformation and exploiting (|44p we may rewrite 
(1451 as 



dt 



N(E-Ei) 



(48) 

where N(E — E m ) is the dimension of the environmen- 
tal band with energy E b = E — E m . This form reflects 
the fact that the dynamics of the occupation probabilities 
with different overall energies are decoupled. We, further- 
more, find from ([48]) for the equilibrium values P®(t — * 
oo) oc N(E — Ei). Thus, the equilibrium state is in accord 
with the a priori postulate in that sense that the probabil- 
ity to find the full system in some subspace is proportional 
to the dimension of this subspace. However, it is in general 
impossible to transform (|4"8|) in a closed set of equations 
for the occupation probabilities pa = J2e Pf 01 ^ alone. 
This may only be done if either only one energy subspace 
E is occupied at all, or if the transition rates ^f m are 
independent of E and the number of states of the envi- 
ronmental bands N a scales as oc exp(/3S a ). Then (|48)l 
may be summed over E yielding 



dt 



Pu(t) 



E 



(*) 



} 7imPii{t)] (49) 



which is the usual closed form for the dynamics of the 
pa with the standard canonical equilibrium state pu(t — * 
oo) oc exp(— (3Ei). Thus it is essentially the exponentially 
growing density of states of typical infinite environments 
that allows for a closed dynamical description of the con- 
sidered system S alone and produces the standard Gibb- 
sian equilibrium state. 



3.5 Thermalization and Decoherence 

In order to investigate the relation between the decay of 
diagonal and off-diagonal elements of the reduced den- 
sity operator of S, we concretely analyze as an example 
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|1> 



AE 



]Se N 3 
JSe N 2 
= J Se Ni 

Fig. 4. Three band model for the investigation of the relation 
between decoherence and thermalization 



V 




Fig. 5. Interaction matrix for the model Fig. [3] canonical 
blocks (gray), microcanonical blocks (white), in case of weak 
coupling irrelevant blocks (hatched). 

a slightly modified model featuring the above mentioned 
structure (exponential state density, equal rates) yield- 
ing autonomous dynamics for p. The model is depicted in 
Fig.Hl For simplicity we consider only three environmental 
bands with the same density of states, i.e., the exponential 
prefactor from (|49|) vanishes, (3 = 0. (This eventually im- 
plies infinite temperature.) As mentioned the rates that 
control the dynamics of the diagonal elements (canoni- 
cal dynamics, thermalization) have to be equal, thus, we 
choose Aoi,i2 = Aoi,23 = A can - The rates that control the 
dynamics of the off-diagonal elements (microcanonical dy- 
namics, decoherence) also have to be equal among them- 
selves, but may differ from the "canonical rates". Thus, 



we choose A, 



00,33 



An 



An 



An 



An 



Since 



all other parts of the interaction would not fulfill the res- 
onance condition anyway, we set them to zero (cf. Fig. [5]) . 
Plugging those model parameters into (|49]) yields 



df 



Me 



(Pll ~ Poo) 



dt Me 
Defining the thermalization time as 



Poo ~ Pn 



T, 



Me 



4ttA2 iV 



(50) 
(51) 

(52) 



the solution of the above set of differential equations is 
just an exponential decay according to e~*/ Tth . 



Apparently the Pjj l0 do not "mix" with respect to 
different a [cf. (146|) ]. Thus, if initially the environment 
only occupies, e.g., band 2, for the full dynamics the off- 
diagonal element of p will be simply given by p+j = fy,2 
[cf. (H2J)]. In this case we find from (|4"rj|) 



djpiol 
dt 



2kN (A^ an + A^ ic 



Me 



\Pw\ 



(53) 



Using the definition £ = A m i C /A can , we find for the deco- 
herence time, i.e., the time-scale on which |pio| decays 



Me 



Me 



2^V(A c \ n + A^ ic ) 2^VA c 2 an (l + £ 2 ) 
2T th 
l + £ 2 ■ 



(54) 



For the absence of microcanonical coupling terms (A m ic = 
— * £ = 0) we get 2T t h = Td cc which is a standard re- 
sult in the context of atomic decay, quantum optics, etc. 
Nevertheless, for increasing £ decoherence may become ar- 
bitrarily faster than thermalization which is a central fea- 
ture of models that are supposed to describe the motion 
of particles subject to heat baths, like, e.g., the Caldeira 
Legget model. Thus our model exhibits a continuous tran- 
sition between those archetypes of behavior. 



4 Application 

4.1 Relaxation Dynamics in Model Systems 

In this Section concrete models are introduced and the 
corresponding time dependent Schrodinger equations are 
solved. Then the results are compared to predictions from 
HAM and standard open system methods. Our first model 
is of the type depicted in Fig. [1] The two level system 
features a splitting of AE = 25w. Here and in the fol- 
lowing we use an arbitrary energy unit u. The environ- 
ment consist of two bands of width Se — 0.5m with the 
same amount of levels N = Ni = N2 = 500 in each one 
and separated also by AE = 25u. As already mentioned 
we use complex Gaussian random matrices to model the 
coupling, thus, satisfying the criteria Sect. 12.11 First, we 
choose only a canonical interaction due to the coupling 
strength A can = 5 • 10~ 4 u (A mic = 0). 

At first we analyze the decay behavior of two differ- 
ent pure product initial states. The environmental part of 
both initial states is a pure state that only occupies the 
lower band, but is apart from that chosen at random. Ir- 
respective of its pureness only with respect to occupation 
numbers, E's initial state can be considered an approxi- 
mation to a Gibbs state with 5e <C &Te <C AE and the 
temperature of the environment Te (in the example at 
hand, e.g., kT^ ~ 5m). For small Se the temperature may 
be arbitrarily small. Initially, the system S is firstly cho- 
sen to be completely in its excited state and, secondly, 
in a 50:50 superposition of ground and excited state. The 
probability [density matrix element pn(t)] to find the sys- 
tem excited as produced by the first initial state is shown 
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Fig. 6. Evolution of the excitation probability of a product 
initial state. Dashed line refers to the standard master equation 
in born approximation. 




1000 2000 
t[H/u] 

Fig. 7. Off-diagonal element for a correlated initial state. 



in Fig. [SI Since the first initial state does not contain any 
off-diagonal elements, we find |poi| 2 ~ for all times. 
This is different for the second initial state investigated in 
Fig. [3 it starts with |poi| 2 = 0-25 and is thus well suited 
to study the decay of the coherence. (The diagonal ele- 
ments of the second state are already at their equilibrium 
value /On(O) = 0.5 in the beginning and exhibits no further 
change.) 

By numerically solving the time-dependent Schrodin- 
ger equation for the full model's pure state we find for the 
reduced state of the system, an exponential decay, up to 
some fluctuations as depicted in Fig. [5] and Fig. [JJ (For the 
baths initial state being a real mixed Gibbs state one can 
even expect fluctuations to be smaller, since fluctuations 
corresponding to various pure addends of the Gibbs state 
will partially cancel each other.) The solid lines are the 
HAM results as computed from (j45|) and (j46|) . Obviously, 
they are in accord with the exact result. 

The full model is Markovian in the sense that bath 
correlations decay much faster than the system relaxes, 
concretely bath correlations decay on a time scale of r c sa 
h/Se = 2 (all times given in units of h/u), whereas the sys- 
tem relaxes on a timescale T\ k, 640 (cf. Fig.[S]). Neverthe- 
less, S's excitation probability deviates significantly from 
what the standard methods (BA) predicts (cf. Fig. [6|): 
The beginning is described correctly, but rather than end- 
ing up at temperature T = Te as the BA predicts for 
thermal environment states Q , S ends up at temperature 
T = oo, i.e., equal occupation probabilities for both levels. 
The equilibrium value of S's excitation probability is given 



ttS, 



400 



200 





Schrodinger x 




HAM — 







1 2 3 4 5 

Fig. 8. Dependence of the decoherence time on the micro- 
canonical interaction strength. HAM theory according to (|54[1 . 



by pn(oo) = JVi/(JVi + N 2 ). Thus, only if N 2 > Ni (infi- 
nite bath) the BA produces correct results. Note, however, 
that it is not the finite density of states that causes the 
break down of the BA, since the BA produces wrong re- 
sults even for Ni , N 2 — > oo as long as the above condition 
is not met. 

Furthermore, a condition often attributed to the BA, 
namely that S and E remain unentangled, is not fulfilled: 
When S has reached equilibrium the full system is in a 
superposition of |S in the excited state (g> E in the lower 
band) and \S in the ground state (g> E in the upper band). 
This is a maximum entangled state with two orthogonal 
addends, one of which features a bath population corre- 
sponding to Te ~ 0, the other a bath population inver- 
sion, i.e., even a negative bath temperature. These find- 
ings contradict the concept of factorizability, neverthe- 
less, HAM predicts the dynamics correctly. This is in ac- 
cord with a result from [§, [ToL ITsl ] claiming that an evo- 
lution towards local equilibrium is always accompanied 
by an increase of system-bath correlations. However, the 
off-diagonal element evolution coincides with the behavior 
predicted by the BA. Thus, in spite of the systems finitc- 
ness and the reversibility of the underlying Schrodinger 
equation S evolves towards maximum local von Neumann 
entropy (see Fig. [5] and Fig. [7]) which supports the con- 
cepts of [12 1 . 

To show that it is indeed possible to get different time 
scales for the decay of diagonal elements of the density 
matrix (thermalization) and the decay of off-diagonal ele- 
ments (decoherence) according to pure Schrodingerian dy- 
namics we consider the concrete model as addressed in 
Sect. 1331 with parameters N = 500, 5e = 0.5u, AE = 25u 



and A r 



5-10- 



However, we choose the microcanon- 



ical interaction strength A m i C , in units of the canonical one 
between £ = and £ = 5. As an initial state we prepared 
a 90:10 superposition of ground and excited state in the 
system, environment somewhere in the middle band. This 
refers to a finite off-diagonal element in the beginning. We 
have computed the Schrodinger dynamics of both diagonal 
and off-diagonal elements of the two level system. By fit- 
ting an exponential to the off-diagonal element we get the 
decoherence time Td cc in dependence of the microcanon- 
ical coupling strength. In Fig. [8] we show this numerical 
decoherence time Td oc in comparison with the theoretical 
prediction of the HAM theory, thus (|54|) . As can be seen, 
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Fig. 9. Deviation of the exact evolution of the spins excitation 
probability from the HAM prediction for a set of entangled 
initial states. 
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Fig. 10. Deviation of the exact evolution of the spins excita- 
tion probability from the HAM prediction for increasing num- 
ber N of states in the environment. 



the numerical result is in very good accordance with our 
theory. 



4.2 Accuracy of HAM 



Since HAM is just a "best guess theory" the exact evo- 
lution follows its predictions with different accuracies for 
different initial states, even if all conditions on the model 
are fulfilled. To analyze this for, say pn(t), we introduce 
D 2 , being the time-averaged quadratic deviation of HAM 
from the exact (Schrodinger) result 

° 2 = ik C dt ( p " AM(<) - ^ act w) 2 • ( 55 ) 

Thus, D is a measure of the deviations from a predicted 
behavior. The results of the investigation for our model 
(Fig. Q]) are condensed in the histogram (Fig. [5J v = 3, 
N = 500). The set of respective initial states is character- 
ized by a probability of 3/4 for |S in its excited state (g> E 
in its lower band) and 1/4 for |S in its ground state (£> E in 
its upper band) . Within these restrictions the initial states 
are uniformly distributed in the corresponding Hilbert 
subspace. Since all of them are correlated the application 
of a product projection operator technique would prac- 
tically be unfeasible. However, as Fig. [9] shows, the vast 
majority of them follows the HAM prediction quite closely, 
although there is a typical fluctuation of D = \/2 ■ 10~ 2 
which is small compared to the features of the predicted 
behavior (which are on the order of one), due to the finite 
size of the environment (cf. also fluctuations in Fig. [6]) . 

In Fig. |TD] the dependence of D 2 on the number of 
states of E is displayed for N = 10, ... , 800 (one evolution 
for each size of the environment). At N — 500 like used in 
the above accuracy investigation we find the same typical 
fluctuation, whereas for smaller environments the typical 
deviation is much bigger. We find that the squared devia- 
tion scales as 1/N with the size of the environment, thus, 
making HAM a reasonably reliable guess for many-state 
environments. 



5 Limits for the Applicability 

The dynamical considerations of Sect. [3] are only guesses, 
but as guesses they are valid for any initial state regard- 
less of whether it is pure, correlated, entangled, etc. Thus, 
in contrast to the standard Nakajima-Zwanzig and TCL 
methods HAM allows for a direct prediction of the typi- 
cal behavior of the system. Nevertheless, for deriving the 
above HAM rate equations we have claimed (and already 
discussed) that there is a reasonably well defined corre- 
lation time t c (cf. Sect. I3.2[) at all. Additionally, we used 
two further approximations: The truncation of the Dyson 
series in second order [see (|22|)] and the replacement of 
the actual value of an expectation value by the average 
in the respective Hilbert space compartment [see (1331) ]. 
In the following we will investigate the validity of these 
approximations in more detail. 



5.1 Truncation of the Dyson Series 

In (|22|) we truncated the Dyson series arguing that for 
short times r and small interaction strength this can be 
a reasonable approximation. We require, however, r > r c . 
Thus, for given interaction strength, the time for which the 
truncation should hold, should exceed the correlation 
time, i.e., > r c . How can be at least approximately 
determined? 

Consider the deviation \6ip(t, t)) of a state at time t+r 
from the state at time t, i.e., \Si/j(t,r)) :— \if)(t + r)) — 
\ip(t)) and let the norm of this deviation be denoted as 
A(t, t) — (5ip(t, T)\5ip(t, r)). If we now evaluate A(t,r) 
by means of a truncated Dyson series and find it small 
compared to one it is consistent to assume that higher 
orders are negligible for the description of \ip(t + r)). If 
we, in contrary, find it to be large compared to one, the 
truncation is definitely not justified. Thus, we implicitly 
define roughly as A(t, to) w 1. 

Truncating the Dyson series to leading order yields (cf. 
A PP . A) 

A{t,T) = m)\U 2 {t,r)\m)- (56) 

Since we in general do not know \ip(t)) in detail, but 
only the P's we replace, following again the argument in 
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Sect. 13. H the actual value of A(t, r) by its Hilbert space to evaluate this for S =: Pij, a {t + r) — Pij. a (t) under the 



average {A(t, ■ 



l {(4>\^J, a \4>)=P*3,a(t)}> 
P, 



thus, obtaining 



A & r ) * E ir Tt ^-^ r)} • (57) 



Exploiting ([86]) we find 



A(t,r) = ^2Re/ lm , Qi) 



i /nab 



(58) 



which, taking (|43|) and (|44|) into account and for times 
r > t c eventually yields 



E 

imab 



-Pii,a^fmi,baT 



(59) 



where the second form refers to the notation introduced 
in and below (|47f . Thus, Z\(i,r) grows linear in r. Since 
all the probabilities Pf(t) sum up to one at all times the 
growth is essentially determined by the rates 7^. Since 
already the sum of the Pf(t) over i and some fixed overall- 
energy is a constant of motion, rates belonging to energy 
subspaces E which are not occupied in the beginning will 
never influence A(t, r). Hence one should consider , the 
time for which the truncation of the Dyson series holds 
within the invariant energy subspace E. From (|59[) we find 
as an rough estimate for 



E 



(60) 



where N$ is the number of eigenstates of S. Comparing 
this to (|45|) and (|46|) it becomes obvious that this is also 
roughly the time-scale for the relaxation dynamics of the 
P's. This implies that the claim t¥ > t c , which guaran- 
tees the applicability of the truncation of the Dyson series, 
is equivalent to claiming that the typical relaxation time 
of the P's should be long compared to the typical corre- 
lation time r c . The latter has already been claimed before 
([7]) in order to transform the iteration scheme into a dif- 
ferential equation. This condition can easily be controlled 
by changing the overall interaction strength A. We find 
that for values of A that violate the above condition the 
agreement between the numerical solution and the HAM 
prediction vanishes. 



5.2 Hilbert Space Variance 

Here we quite briefly consider the assumption that gave 
raise to the replacement of actual expectation values by 
their Hilbert space averages in Sect. 13.11 As already men- 
tioned, such a replacement can only yield a reasonable 
result if the largest part of the possible expectation values 
is indeed close to the corresponding Hilbert space average. 
To analyze this we consider the Hilbert space variance of, 
aayA^\SW^.e.,A H S=MS\^) 2 j-msm 2 .HA H S 
is small the above condition is satisfied. We would like 



restriction of given P^ ja (i). This, however, turns out to 
be mathematically rather involved and we have not man- 
aged to do so, yet. But, for the Hilbert space variance of 
any Hermitian operator S without any restriction one gets 
(cf. 0) 



A H S 



1 



N + l 



Tr{S 2 } /Tr{S} 



N 



N 



(61) 



where the term in brackets obviously is the spectral vari- 
ance of S and N denotes the dimension of the full sys- 
tem. At this point it simply appears plausible (which is 
of course far from being a proof) that the spectral vari- 
ances of the above defined S remain constant if one varies 
N, but keeps the rates 7 constant. Thus, for growing TV 
the replacement becomes more and more justified. Such 
a scenario is in accord with the general ideas of quantum 
thermodynamics as presented in [15| and especially backed 
up by the numerical findings of Sect. 14.21 



6 Conclusion 

Explicitly exploiting the Hilbert Space Average Method 
(HAM) we have in essence shown, that statistical relax- 
ation may emerge directly from the Schrbdinger equation. 
This requires the respective system being coupled in an 
adequate way to a suitable environment. This environ- 
ment must feature many eigenstates. There is, however, no 
minimum particle number limit. Thus the thermodynamic 
limit appears to be essentially controlled by the number of 
environmental eigenstates involved in the dynamics rather 
than by the number of environmental particles. This re- 
laxation behavior results even for correlated initial states, 
nevertheless, standard open system methods may fail to 
produce the correct result. 



We are indebted to H.-P. Breuer and G. Mahler for interesting 
discussions on this subject. Financial Support by the Deutsche 
Forschungsgemeinschaft is gratefully acknowledged. 



Appendix A: Time evolution 



As already mentioned in Sect. 12.11 in the interaction pic- 
ture, V itself earns a time dependence 



V(t) = QiHioct/fr y e -iH loc t/h _ 



(62) 



However, due to the organization of the interaction ([5]) 
the total time dependence may be assigned only to the 
environment parts 



(63) 



with 



C l3 {t) = e^'/^e-^e^-Wi . (64) 
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[We already mention this here for later reference, cf. (|75|) ] . 
Assuming weak interactions ()17|) and short times r the 
time evolution D(r, t) resulting from the Dyson series may 
be truncated at second order, i.e., 

m + r)) « [l - -U!^) - ^U 2 (T,t)]m)) > (65) 

V v ' 

D 2 (r,t) 

with the two time evolution operators 



66|) and the interaction (|63|) one gets 



Ui(r,t) 



dr'V (r' + t) , 



(66) 



U 2 (t, t) = f T+t dr' f T +t dr" V(t' + t) V(t" + t) 



(67) 



Note that the integration in |67|) is time ordered, i.e., t > 
t' > t" . Furthermore, the first order operator U\{t) is 
Hcrmitian due to the Hermiticity of the interaction. 



Appendix B: Correlation Functions 

One has to analyze the traces on the right hand side (|34|) 
Let us therefore abbreviate those term by S. Using ((65 
we get 

3 = Tr{DlP ijt a D 2 P ifj >, a >} = 

i ^ 1 



Tr 



where we used the Hermiticity of the operator U\. 

To evaluate this complicated trace expression we will 
consider each order of time evolution operators in (|68p 
separately, defining 

S = So + Si + S 2 ■ (69) 
Using (fTTj) the zeroth order of (|68|) yields 

S = Tl-{Pij, a Pi'j',a'} = Si'jSfiSa'aNa . (70) 

By a cyclic rotation within the trace the first order 
may be written as 



"5*1 — ~rT^{UiPij, a Pi' j' a ' Pi' j\a' Pij,aU\\ 

= ^a'aO^T^lAi'.J - 8j'iR{UiPi>j, a }) , (71) 

where we used (|11[) again. Concentrating on the first term, 
introducing the definition of the time evolution operator 



T+t 



TiiUxPij^a} = JdT'TY{V(T' + t)P ir n a } 

fT + t 

= dr' J2 Tr{PkiC k i{T' + t)P lf fl a } 

= fdr' Y^TrsipMPi^TrviCkiiT' +t)n a } 
Jt 



T + t 



dT'Tr E {C fi (T' + t)n a }. 



(72) 



Due to the condition on the interaction @ those terms 
are zero. We find an analogous result for the second trace 
of (|71jl and thus we finally end up with 

Si = . (73) 

For the second order terms of (|68|) we get 

S 2 =Jp^{UlPij,a,U\Pitji, a i 

- Sji'da'aU^Pij'.a' ~ Sj>i5a>aU 2 Pi>j,a>} ■ (74) 

We concentrate first on the last term, plugging in the def- 
inition of U 2 from ([FT]) yields 

rT + t p T '+t 

Tr{U2Pi'j,a>} = J t dT J t dT " x 

x Tr{U(r' + t)V{r" + t)A',,.'}. ( 75 ) 

exploiting (|63|) and performing the trace with respect to 
S we find 

= fdr' fdr" Tr E {Q m (r' + t)C mi > (r" + t)fl a ,} . 
Jt Jt 

(76) 

Since the operators that generate the time-dependence of 
V(t) [cf. (f6"4")) ] commute with TI a i and due to the invari- 
ance of the trace with respect to cyclic permutations of the 
traced operators, the above "projected correlation func- 
tions" only depend on the difference between the time ar- 
guments of the U's. Since then the integrand no longer de- 
pends on t, the t which appears in the integration bound- 
aries may simply be set to zero. Hence one finds for the 
above expression 

= V/ dr' / dT"Tr E {C jm (T' -T")C m ,n a ,}. (77) 
m Ja Jo 

As argued in the beginning the parts of the interaction are 
uncorrelated unless they are not adjoints of each other. 
This means that the above traces can only be nonzero for 
the case j — i' . Furthermore, one does the transformation 
(r' - t") -> r", thus, 



, ; / dr' / dr"Tr E {C jm {T")C mj n a/ }. (78) 
o Jo 
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Finally, plugging in the unit operator of the environment 
in terms of projection operators, one finds 

J o Jo b 
= Y,k'j I dr ' / d T"Tr E {C jm , a , b (T")C mjM ,} . 

rnb J ° J ° 

(79) 

Comparing this to (137|) we end up with 

Tr{U2Pi>j,a>} = 6i 'i fomM r ) ■ ( 8 °) 



same integral transformation described before (|78p . even- 
tually write 

T±{UiPii,aUiPi'j',a'} = ^'2Re /«»,„. (r) . (86) 

Putting all the bits and peaces from l[70]). (|73| . (|74|) . 
(|50). (|5T]) and (dHD together, we eventually find 

Tr{DlPij, a D 2 Pi>j> : a>} = Si'jSfiSa'aNa 

Sij Si'ji 2Re/ii/ !0a /(r) 

" 5 i fi Sa ' a [f*™,a>b( T ) + fjm,a'b(T)] ) • (87) 



Completely analogous we find for 

TrtflPij',*'} = £ S i'i f!m,a'b(r) ■ (81) 



rnb 



It remains the computation of the first term of (|74|) . 
Using the same argumentation as before (the fact that 
there are no correlations between different parts of the 
interaction as well as a cyclic rotation within the trace 
operation) we find for the trace 

/■T+t rT+t 

Jt Jo 



X 



Tr E {C u >(T" + t)n a ,C Vi {r' + *)#«} . (82) 



By the same arguments which are given below (|75[) this 
may be written independently of the absolute time t 

= SiiSi, r I dr' f dT"Ti^{C lt> {T" -T')fl a ,C tli na\. 



JO 



(83) 



The (non-time-ordered) integration of the above expres- 
sion may be written in terms of a time-ordered integration 
by adding the time-reversed integrand 



dr' f dr"Tr E {CW(r" - r')/7 'C«^ } 
o Jo 

+ f dr' / dr"Ti E {Cw(T' -T")n a/ C^n a }. 
Jo Jo 



(84) 



Shifting in the second term the time dependence to the 
other C operator and expressing everything within the 
trace by its adjoint yields 



T PT 

I i II r 



= dr' dT"Ti E {C u i(T " -r')U a ,C Vi n a } 
Jo Jo 

+ f dr' / dr"Ti E {(Cu'(r" -T^n^C^flay}- (85) 
Jo Jo 
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